It is a testrun for v3b model checking recovery.
v3b1 model was used with standard setup (slightly wider prior on sigma + stan-specific phi calculation)
Runtime ranged 2.1-2.7k sec
3 divergences
pairs( fit, pars = c( paste0( 'r_lower[', 1:10, ']'), 'lp__'))
pairs( fit, pars = c( paste0( 'r_upper[', 21:30, ']'), 'lp__'))
pairs( fit, pars = c( 'sigma_lower', 'sigma_upper', 'lp__'))
A few divergences and funnel still clearly visible, but on the overall, it appears acceptable
pairs( fit, pars = c( 'sigma_gains', 'lp__'))
pairs( fit, pars = c( 'sigma_lpa', 'lp__'))
pairs( fit, pars = c( 'sigma_losses', 'lp__'))
pars <- c( 'sigma_lower', 'sigma_upper')
mcmc_recover_intervals( posterior[,,pars], true_vals[true_ix,pars])
Kind of overestimated standard deviations here
pairs( fit, pars = c( 'theta_upper', 'lp__'))
pairs( fit, pars = c( 'theta_lower', 'lp__'))
mcmc_parcoord( posterior[,,c( 'sigma_lower', 'sigma_upper')], np = np)
mcmc_parcoord( posterior[,,grep( 'sigma_gains', dimnames( posterior)[[3]])], np = np)
mcmc_parcoord( posterior[,,grep( 'sigma_lpa', dimnames( posterior)[[3]])], np = np)
mcmc_parcoord( posterior[,,grep( 'sigma_gains', dimnames( posterior)[[3]])], np = np)
mcmc_parcoord( posterior[,,grep( 'L_lpa', dimnames( posterior)[[3]])], np = np)
mcmc_parcoord( posterior[,,grep( 'L_gains', dimnames( posterior)[[3]])], np = np)
mcmc_parcoord( posterior[,,grep( 'L_losses', dimnames( posterior)[[3]])], np = np)
lapply( c( 'lpa', 'gains', 'losses'), function( dd) {
apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
c( 1, 2), function( x) {
x %>%
matrix( ncol = sqrt( length( .))) %>%
`%*%`( t( .)) %>% `[`( lower.tri( .))
}) %>% aperm( perm = c( 2, 3, 1)) %>%
`dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
mcmc_parcoord( np = np)
})
## [[1]]
##
## [[2]]
##
## [[3]]
sapply( c( 'gains', 'lpa', 'losses'),
function( dd) {
apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
c( 1, 2), function( x) {
x %>%
matrix( ncol = sqrt( length( .))) %>%
`%*%`( t( .)) %>% `[`( lower.tri( .)) %>% median()
})
}, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
mcmc_parcoord( np = np)
sapply( c( 'gains', 'lpa', 'losses'),
function( dd) {
apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
c( 1, 2), function( x) {
x %>%
matrix( ncol = sqrt( length( .))) %>%
`%*%`( t( .)) %>% `[`( lower.tri( .)) %>% abs() %>% median()
})
}, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
mcmc_parcoord( np = np)
sapply( c( 'gains', 'lpa', 'losses'),
function( dd) {
apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
c( 1, 2), function( x) {
x %>%
matrix( ncol = sqrt( length( .))) %>%
`%*%`( t( .)) %>% `[`( lower.tri( .)) %>% abs() %>% log() %>% sum()
})
}, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
mcmc_parcoord( np = np)
ff <- function( x) { sum( log( abs( x))) * prod( sign( x))}
sapply( c( 'gains', 'lpa', 'losses'),
function( dd) {
apply( posterior[,,grep( paste0( 'L_', dd), dimnames( posterior)[[3]])],
c( 1, 2), function( x) {
x %>%
matrix( ncol = sqrt( length( .))) %>%
`%*%`( t( .)) %>% `[`( lower.tri( .)) %>% ff()
})
}, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], parameters = list( paste0( 'omega_', 1:( dim( .)[3]))))) %>%
mcmc_parcoord( np = np)
rm( ff)
sapply( c( 'u', 'l'), function( dd) {
posterior[,,grep( paste0( '^r_', dd), dimnames( posterior)[[3]])] %>%
apply( c( 1, 2), function( x) median( abs( x)))
}, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], list( parameters = c( 'upper', 'lower')))) %>%
mcmc_parcoord( np = np)
sapply( c( 'u', 'l'), function( dd) {
posterior[,,grep( paste0( '^r_', dd), dimnames( posterior)[[3]])] %>%
apply( c( 1, 2), function( x) mean( abs( x)))
}, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], list( parameters = c( 'upper', 'lower')))) %>%
mcmc_parcoord( np = np)
sapply( c( 'u', 'l'), function( dd) {
posterior[,,grep( paste0( '^r_', dd), dimnames( posterior)[[3]])] %>%
apply( c( 1, 2), function( x) min( abs( x)))
}, simplify = 'array') %>%
`dimnames<-`( c( dimnames( .)[1:2], list( parameters = c( 'upper', 'lower')))) %>%
mcmc_parcoord( np = np)
pars <- dimnames( posterior)[[3]]
pars <- pars[!grepl( 'pe_lpa', pars)]
mapply(function( x, a, b) (x > a) & (x < b),
x = true_vals[true_ix,pars],
a = summary( fit)$summary[pars,'2.5%'],
b = summary( fit)$summary[pars,'97.5%']) %>% mean()
## [1] 0.9193254
Somewhat lower, but still acceptable recovery of parameter values (91%)